#include<../data_tools/usefull_routines.c>

int main(int argc, char *argv[])
{
  int i,j,ndisk,nbines,npart,nreal=0,status;
  float zmin,zmax,deltaz,mass,zini,zfin,line; 
  float meanmass,*z,*sigma,*errors,f,p;
  double *Z,size;
  char *infile,outfiles[100];
  size_t *ind;
  particulas *diskpart;
  infile = argv[1]; 
  nbines = atoi(argv[2]); 

  printf("\n");
  printf("#############################################################\n");
  printf("Computing vertical density profile for the disk:\n");
  printf("%s\n",argv[1]);
  printf("#############################################################\n");
  printf("\n");

  //Reading datas
  FILE *diskfile;
  diskfile = fopen(infile,"r");
  ndisk = contadorParticulas(infile);
  printf("Reading %d particles\n",ndisk);
  
  diskpart = (particulas *)malloc((size_t)ndisk*sizeof(particulas));
  if(diskpart == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }

  Z = (double *)malloc((size_t)ndisk*sizeof(double));
  if(Z == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }

  z = (float *)malloc((size_t)nbines*sizeof(float));
  if(z == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }
  sigma = (float *)malloc((size_t)nbines*sizeof(float));
  if(sigma == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }
  errors = (float *)malloc((size_t)nbines*sizeof(float));
  if(errors == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }

  ind= (size_t *)malloc((size_t)ndisk*sizeof(size_t));
  if(ind == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }

  for(i=0; i<ndisk; i++)
    {
      fscanf(diskfile,"%d %e %e %e %e %e %e %e %e\n",
	     &diskpart[i].id,
	     &diskpart[i].pos[0],&diskpart[i].pos[1],&diskpart[i].pos[2],
	     &diskpart[i].vel[0],&diskpart[i].vel[1],&diskpart[i].vel[2],
	     &diskpart[i].mass,
	     &diskpart[i].pot);
      Z[i] = diskpart[i].pos[2];
      ind[i] = i;
    }

  //Sorting Z anf finding zmax 
  gsl_sort_index(ind,Z,1,ndisk);
  
  if(fabs(Z[ind[ndisk-1]]) >= fabs(Z[ind[0]]));
  {
    zmax = fabs(Z[ind[ndisk-1]]);
    zmin = -1.0*zmax;
  }
  if(fabs(Z[ind[ndisk-1]]) < fabs(Z[ind[0]]));
  {
    zmax = fabs(Z[ind[0]]);
    zmin = -1.0*zmax;
    }
  printf("Done.\n");
  printf("\n");

  zmax =6.0;
  zmin = -zmax;

  //Computing vertical density profile
  printf("Computing vertical density profile\n");
  FILE *perfilfile;
  sprintf(outfiles,"%s_vertical_density.dat",infile);
  perfilfile = fopen(outfiles,"w");

  deltaz = (zmax-zmin)/nbines;
  zini = zmin;

  for(i=0; i<nbines; i++)
    {
      npart = 0;
      mass = 0.0;
      zfin = zini + deltaz;

      for(j=0; j<ndisk; j++)
	{
	  if((diskpart[j].pos[2]>=zini) && (diskpart[j].pos[2]<zfin))
	    {
	      npart = npart + 1;
	      mass = mass + diskpart[j].mass;
	    }
	}
      line = zfin-zini;
      meanmass = mass/npart;
      z[i] = (zini+zfin)/2.0;
      sigma[i] = mass/line;
      errors[i] = meanmass*sqrt(npart)/line;
      if(npart==0) 
	{
	  z[i] = 0.0;
	  sigma[i] = 0.0;
	  errors[i] = 1.0;
	  nreal = nbines-1;;
	}

      if(npart>5)
	{
	  fprintf(perfilfile,"%e %e %e\n",z[i],sigma[i],errors[i]);
	}
      zini = zfin;
    }
  printf("Done.\n");
  printf("\n");
 
  //Fitting vertical density profile
  printf("Fitting vertical density profile\n");
  double densidadVertical(double z,double A,double zo,double center)
  {
    double Cosh,perfil;
  
    Cosh = 1.0/(cosh((z-center)/zo));
    
    perfil = A*(1.0/zo)*(pow(Cosh,2));
    
    return perfil;
  }
  //funcion que calcula el chisquare
  double ChiSquare(const gsl_vector *parmodelo, void *param)
  {
    int i;
    double a,center,zo;
    double chisquare;
    double xo,yo,s,yt;
    
    a = gsl_vector_get(parmodelo,0);
    zo = gsl_vector_get(parmodelo,1);
    center = gsl_vector_get(parmodelo,2);
    
    //CALCULO DEL CHI-SQUARE
    chisquare = 0.0;
    for(i=0; i<nbines; i++)
      {
	xo=z[i];
	yo=sigma[i];
	s=errors[i];
	yt=densidadVertical(xo,a,zo,center);
	//Calculo del chisquare
	chisquare+=gsl_pow_int(yt-yo,2)/(s*s);
      }
    
    return chisquare;
  }

  //MINIMIZACION DEL CHISQUARE
  gsl_vector *parmodelo=gsl_vector_alloc(3);
  gsl_vector *dpar=gsl_vector_alloc(3);
  
  gsl_vector_set(parmodelo,0,1.0);
  gsl_vector_set(parmodelo,1,1.0);
  gsl_vector_set(parmodelo,2,0.0);
   
  gsl_vector_set_all(dpar,0.1);
  
  //1-DECLARACION
  const gsl_multimin_fminimizer_type *tipo;
  gsl_multimin_fminimizer *minimizador;
  gsl_multimin_function funcion;
  
  //2-ALOCACION E INICIALIZACION
  tipo=gsl_multimin_fminimizer_nmsimplex;
  minimizador=
    gsl_multimin_fminimizer_alloc(tipo,3);
  
  funcion.f=&ChiSquare;
  funcion.n=3;
  funcion.params=NULL;
  
  gsl_multimin_fminimizer_set(minimizador,
			      &funcion,
			      parmodelo,
			      dpar
			      );
  
  //3-USO
  do{
    gsl_multimin_fminimizer_iterate(minimizador);
    //Es un estimativo de la distancia al minimo
    size=gsl_multimin_fminimizer_size(minimizador);
    status=gsl_multimin_test_size(size,1E-5);
  }while(status!=GSL_SUCCESS);
  printf("Done.\n");
  printf("\n");  

  parmodelo=gsl_multimin_fminimizer_x(minimizador);
    fprintf(stdout,"zo = %e\n",
	  gsl_vector_get(parmodelo,1));
  fprintf(stdout,"zmax = %e\n",
	  gsl_vector_get(parmodelo,2));
 fprintf(stdout,"coef = %e\n",
	  gsl_vector_get(parmodelo,0));
 fprintf(stdout,"Chisquare = %e with %d points\n",
	 gsl_multimin_fminimizer_minimum(minimizador),nreal);
 printf("Done.\n");
 printf("\n");

  FILE *fitfile;
  sprintf(outfiles,"%s_vertical_density.fit",infile);
  fitfile = fopen(outfiles,"w");
  for(f=zmin; f<zmax; f+=0.01)
    {
      p = densidadVertical(f,gsl_vector_get(parmodelo,0),gsl_vector_get(parmodelo,1),gsl_vector_get(parmodelo,2));
      fprintf(fitfile,"%e %e\n",f,p);
    }

  FILE *scale;
  sprintf(outfiles,"%s.zo",infile);
  scale = fopen(outfiles,"w");
  fprintf(scale,"%e",gsl_vector_get(parmodelo,0));
  fclose(scale);

  //4-LIBERACION
  gsl_multimin_fminimizer_free(minimizador);


  free(diskpart);
  free(Z);
  free(z);
  free(sigma);
  free(errors);
  free(ind);
  
  fclose(diskfile);
  fclose(perfilfile);
  fclose(fitfile);
 
  return 0;
}
